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PREFACE 


Under  the  sponsorship  of  the  Tactical  Technology  Office  of  the 
Defense  Advanced  Research  Projects  Agency,  The  Rand  Corporation  has 
been  engaged  in  the  analysis  and  development  of  hydrodynamic  design 
criteria  employing  concepts  of  boundary-layer  control,  including  shap- 
ing, suction,  and  heating. 

This  report  supplies  a theoretical  explanation  of  flow  develop- 
ment in  a heated  pipe.  The  analysis  shows  that  a secondary  flow  is 
generated  in  the  core  of  the  pipe  by  heating  the  wall,  and  that  this 
flow  differs  from  that  over  the  exterior  of  bodies  of  revolution.  The 
results  define  domains  where  this  secondary  flow  can  be  neglected  and 
where  test  data  can  be  used  in  interpreting  exterior  flow.  They  are 
applied  to  a tube  test  set  up  at  Colorado  State  University,  Fort 
Collins,  Colorado,  by  the  Autonetics  Division  of  Rockwell  Interna- 
tional, to  experimentally  verify  theories  predicting  stabilization  of 
the  boundary  layer  by  heating,  shaping,  and  other  means. 

This  report  should  be  useful  to  hydrodynamlcists , to  designers 
of  submerslbles , and  to  others  interested  in  applying  fluid  mechanics 
to  improve  the  performance  of  underwater  vehicles . 

Other  related  Rand  publications  include: 

R-1752-ARPA/ONR,  Low-Speed  Boundary -Layer  Transition  Workshop, 

W.  S.  King,  June  1975. 

R-1789-ARPA,  Controlling  the  Separation  of  Laminar  Boundary 
Layers  in  Water:  Heating  and  Suotion,  J.  Aroesty  and  S.  A. 

Berger,  September  1975. 

R-1863-ARPA,  The  Effects  of  Wall  Temperature  and  Suction  on 
Laminar  Boundary -Layer  Stability,  W.  S.  King,  April  1976. 

R-1898-ARPA,  "e^”:  Stability  Theory  and  Boundary -Layer  Transi- 

tion, S.  A.  Berger  and  J.  Aroesty,  February  1977. 

R-1907-ARPA,  Buoyancy  Cross-Flow  Effects  on  the  Boundary  Layer 
of  a Heated  Horizontal  Cylinder,  L.  S.  Yao  and  I.  Catton, 

April  1976. 

R-1966-ARPA,  The  Buoyancy  and  Variable  Viscosity  Effects  on  a 
Water  Laminar  Boundary  Layer  Along  a Heated  Longitudinal 
Horizontal  Cylinder,  L.  S.  Yao  and  I.  Catton,  February  1977. 
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SUMMARY 


The  developing  flow  in  the  entry  region  of  a heated  horizontal 
pipe  is  analyzed.  Due  to  the  similarity  of  flow  phenomena  between 
heated  straight  pipes  and  unheated  curved  pipes,  the  techniques  de- 
veloped in  treating  the  flow  problem  in  the  latter  can  be  applied  to 
study  the  flow  development  in  the  former.  The  asymptotic  solution 
of  the  developing  flow  near  the  entrance  of  the  heated  straight  pipe, 
distance  0(a),  is  obtained  by  perturbing  the  solution  of  the  develop- 
ing flow  in  an  unheated  straight  pipe.  Two  vortices  result  from  the 
combination  of  the  radial-directional  and  the  downward  motions  of  the 
fluid  particles  which  are  Induced  by  the  displacement  of  the  boundary 
layer  and  develop  along  the  pipe.  The  axial  velocity  has  a concave 
profile  with  its  maximum  off  the  center  line  near  the  entrance,  and  it 
grows  toward  a uniformly  distributed  profile  downstream.  The  downward 
stream  caused  by  the  displacement  of  the  secondary  boundary  layer 
forces  the  axial  velocity  profile  to  turn  counterclockwise  continu- 
ously along  the  pipe  if  the  flow  is  from  left  to  right.  A favorable 
pressure  gradient  is  generated  on  the  bottom  wall  of  the  pipe;  an  un- 
favorable pressure  gradient  is  Induced  on  the  top  wall. 

The  results  of  the  analysis  are  applied  to  the  "Colorado  tube 
test"  set.  It  has  been  found  that  under  test  conditions,  when  the 
flow  speed  is  not  fast  enough,  strong  wall  heating  can  cause  the 
boundary-layer  flow  to  separate  from  the  top  wall  of  the  pipe,  caus- 
ing transition  to  turbulence;  in  this  case  the  pipe  flow  cannot  be 
used  to  simulate  an  external  boundary-layer  flow. 
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SYMBOLS 


a = radius  of  the  pipe 

A = variable  viscosity  function,  Eq . (13) 

D = Dean  number , = a ‘Re 
fp,  F^,  ^2  “ dimensionless  stream  functions 

2 

g = gravitational  acceleration  = 32.2  ft/sec 

G = dimensionless  temperature  function 

3 

Gr  = Grashof  number  = 

Iq,  = Bessel  functions 

k = thermal  conductivity 

Nq , N = viscosity  ratio,  Eq.  (13) 

P = pressure 

Pr  = Prandtl  number 

q = wall  heat  flux,  Eq.  (1) 
w 

r,  (j),  z = coordinates.  Fig.  1 (ft) 

r = coordinate  normal  to  the  wall 

Ra  = Rayleigh  number 

Re  = Reynolds  number 

T = temperature  (°F) 

AT  = T - T.  (°F) 
w xn 

T = wall  temperature  (°F) 
w 

T.  = inlet  temperature  (°F) 
xn 

u,  v,  w = velocities  in  the  boundary  layer  (ft/sec) 

U,  V,  W = velocities  of  the  invlscid  core  flow  (ft/sec) 


X,  Y,  z = coordinates.  Fig.  1 
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a = curvature  ratio 
a = constant  viscosity,  Eq.  (13) 
3^,  ^3  “ constants.  Table  1 

y = expansion  coefficient 

6 = 

/Re 

e = Eq.  (7) 

Re 


Subscripts 

0,  01,  10,  11  indicate  the  order  of  the  perturbing  functions 
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I.  INTRODUCTION 


The  laminar  flow  convective  heat  transfer  in  a pipe  has  been  a 
research  topic  for  about  half  a century.  Continuous  efforts  have  been 
devoted  to  the  topic  due  to  its  practical  applications  in  various 
engineering  systems.  It  has  been  found  that  natural  convection  and 
variable  material  properties  play  important  roles  in  the  heat  transfer 
and  the  fluid  flow  in  a heated  pipe.  Frequently  the  prediction  of  the 
heat  transfer  by  forced  convection  alone,  without  considering  the 
secondary  flow,  can  cause  large  errors.  Comprehensive  reviews  of  pre- 
vious work  can  be  found  in  two  recent  papers  by  Siegwarth  et  al.  (1969) 
and  Hong  and  Bergles  (1976).  Of  the  early  works,  three  papers  worthy 
of  special  attention  are  summarized  below. 

Morton  (1959)  found  that  the  flow  in  a heated  straight  pipe  is 
similar  to  that  in  a curved  pipe.  He  followed  the  similar  expansion 
technique  developed  by  Dean  (1927,  1928)  in  solving  the  problem  of  a 
fully  developed  flow  in  curved  pipes  and  was  able  to  give  the  solution 
of  laminar  convection  in  uniformly  heated  horizontal  pipes.  He  showed 
that  the  expansion  parameter  is  the  Rayleigh  number,  Ra  (or  the  Grashof 
number,  Gr) . However,  for  a constant  wall  heat  flux  condition,  the 
heat  transfer  is  enhanced  by  the  speed  of  the  flow  and  the  length  of 
a heated  pipe.  Therefore,  the  flow  velocities  and  the  heat  transfer 
also  depend  on  the  Reynolds  number.  Re.  Consequently,  the  solutions 
of  constant  wall  heat  flux  depends  on  the  product  of  the  Rayleigh  and 
the  Reynolds  numbers,  RaRe.  Since  the  analysis  given  by  Dean  is  limited 
to  the  case  of  small  Dean  number,  D,  the  ratio  of  the  centrifugal  force 
and  the  viscous  force,  the  Morton  solution  is  also  valid  only  when  ReRa, 
the  ratio  of  the  buoyancy  force  and  the  viscous  force,  is  small  in  a 
heated  pipe.  For  the  curved  pipe  flow  of  large  D,  Baura  (1963)  ob- 
served that  viscous  forces  are  important  only  in  a thin  boundary  layer 
near  the  wall  and  that  the  motion  outside  the  boundary  layer  is  mostly 
confined  to  planes  parallel  to  the  plane  of  symmetry  of  the  pipe.  With 
these  observations  and  assumptions,  Baura  obtained  an  asymptotic 
boundary-layer  solution;  the  resistance  coefficient  obtained  from  this 
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theory  agrees  with  the  observations  of  White  and  Adler  (see  Baura,  1963) . 
Mori  and  Futagami  (1967)  presented  a theoretical  model  similar  to  Baura' s 
ideal;  the  Nusselt  number  obtained  from  their  model  agrees  well  with 
their  experimental  data.  Another  interesting  character  of  the  flow  in 
heated  pipes,  revealed  by  Mori  et  al.  (1965)  is  that  the  secondary  flow 
generated  by  heating  can  suppress  the  turbulence  level  when  the  inlet 
turbulence  level  is  high  and  can  enhance  it  when  its  level  is  low  at  the 
inlet.  They  also  found  that  the  critical  Reynolds  number  of  the  flow  in 
heated  pipes  converges  to  a value  approximately  equal  to  3.8  x 10  at 
ReRa  = 5 x 10^  independently  on  the  inlet  turbulence  level.  The  similar 
phenomenon  was  observed  by  Taylor  (1929)  for  the  flow  in  curved  pipes. 

In  this  report  we  study  the  entry-developing  flow  in  a heated  straight 

pipe  under  the  conditions  of  constant  wall  temperature  and  constant  wall 

heat  flux,  respectively.  The  inlet  velocity,  is  uniformly  distributed. 

The  wall  of  the  pipe  is  maintained  either  at  a higher  temperature,  T , than 

w 

the  inlet  fluid  temperature,  or  at  the  constant  heat  flux,  q^.  The 

governing  parameters  of  the  problem  are  (see  Eqs.  (2)  in  Section  II): 


(i) 


(ii) 


(iii) 

(iv) 


the  Reynolds  number 


the  Grashof  number 


the  Prandtl  number 
the  ratio 


Re  = 


W.  a 
in 


Yga  (T^  - T.^) 

Gr  = (constant  wall  temperature) 

V 


4 

Yg^ 

kv^ 


(constant  wall  heat  flux) 


(1) 


Pr  = 


y_ 

K 

Gr 

Re^ 

Gr 


Re 


5/2 


(constant  wall  temperature) 


(constant  wall  heat  flux) 


where  a is  the  radius  of  the  pipe,  v is  the  kinematic  viscosity,  6 is  the 
thermal  expansion  coefficient,  g is  the  gravitational  acceleration,  and 
K is  the  thermal  diffusivity. 
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Since  the  development  of  the  entry  flow  in  heated  pipes  is  similar 
to  that  of  curved  pipes  and  the  solution  we  will  present  is  obtained  by 
perturbing  the  entry  flow  in  unheated  straight  pipes,  we  shall  briefly 
review  previous  work  on  the  entry  flow  in  unheated  straight  pipes  and 
curved  pipes. 


UKHEATED  STRAIGHT  PIPES 

Previous  work  on  developing  flow  in  unheated  straight  pipes  and 
two-dimensional  channels  falls  into  the  following  four  categories: 

(1)  linearization  of  the  momentum  equations;  (il)  two-zone  models,  in 
which  a boundary- layer  flow  matches  with  the  downstream  perturbation 
solution  for  the  fully  developed  flow;  (ill)  momentum- Integral  tech- 
niques; and  (iv)  finite-difference  solutions.  A detailed  review  of 
the  published  papers  falling  into  the  above  four  categories  has  been 
summarized  by  Yao  and  Berger  (1975) . From  a critical  analysis  of  the 
entry  flow  problem.  Van  Dyke  (1970)  pointed  out  that  there  are  two 
length  scales  in  channel  entry  flow:  a,  the  half-width  of  the  channel, 

and  aRe,  Most  of  the  early  work  on  this  problem  is  valid  only  for  down- 
stream distance  O(aRe).  Van  Dyke  obtained  a solution  in  the  upstream 
region,  for  distance  0(a),  and  so  reconciled  the  discrepancy  between  the 
earlier  theoretical  results  and  the  experimental  data  near  the  entrance 
of  the  channel.  Independently,  a similar  but  more  detailed  mathematical 

analysis  was  given  by  Wilson  (1971) . 

CURVED  PIPE 

Two  parameters  govern  the  developing  entry  flow  in  curved  pipes: 
a,  the  curvature  ratio,  and  D = Ae,  the  Dean  number.  The  entry  flow 
in  curved  pipes  can  be  categorized  into  the  following  three  cases: 

1.  When  both  a and  D are  small,  the  flow  is  slightly  distorted 
from  the  entry  flow  in  unheated  straight  pipes,  and  the  solution  can 
be  obtained  by  perturbing  the  solution  of  the  developing  flow  in  unheated 
straight  pipes.  There  are  two  regions:  0(a)  and  O(aRe).  The  solutxon 

in  the  region,  for  distance  0(a),  was  given  by  Singh  (1974).  As  the 
fluid  enters  the  pipe,  a boundary  layer  like  that  in  an  unheated  straight 
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pipe  develops  on  the  wall  but  with  a small  azimuthal  component  of 
velocity  due  to  the  inwardly  directed  pressure  gradient,  the  latter 
arising  from  the  circular  nature  of  the  main  central  flow.  The  dis- 
placement effect  of  the  boundary  layer  in  turn  causes  an  accelera- 
tion of  the  flow  in  the  central  core  plus  a secondary  flow  in  the 
cross-sectional  plane.  Initially  there  is  an  inward  flow  from  the 
entire  pipe  circumference,  with  two  singularities  in  the  central  region, 
a node-like  sink  at  the  origin,  and  a saddle-point-like  stagnation 
point;  the  latter  singularity  moves  outwards  as  the  fluid  moves  down- 
stream, finally  vanishing  as  the  cross  flow  (from  the  inside  to  the 
outside  of  the  bend)  sets  in.  (For  further  details,  plus  the  inter- 
esting consequences  this  flow  picture  has  on  the  location  of  the 
point  of  maximum  shear,  see  Singh  (1974).)  The  flow  development  de- 
scribed above  occurs  within  a distance  0(a)  from  the  entrance.  The 
matched  asymptotic  solution  for  this  region  (Singh,  1974)  breaks  down 
at  a distance  0(a//a),  corresponding  physically  to  the  point  beyond 
which  the  effect  of  the  centrifugal  force,  initially  small,  becomes 
as  important  as  inertia  and  viscous  forces  within  the  boundary  layer. 
Also,  the  matched  asymptotic  solution  developed  for  the  region,  dis- 
tance 0(a),  apparently  breaks  down  at  a distance  O(aRe),  since  aRe  is 
the  length  scale  of  the  region  that  flow  transits  to  the  fully  de- 
veloped flow  in  unheated  straight  pipes.  There  are  two  possible  length 
scales  for  the  region  downstream  of  the  region  distance  0(a)  from  the 
entrance.  However,  aRe  < a/r^  when  D is  small.  This  suggests  that  the 
flow  will  be  fully  developed  before  centrifugal  forces  lose  their 
secondary  role,  when  D is  small. 

2.  For  the  case  when  Ot  is  small  and  D is  large,  the  flow  appears 
to  develop  as  follows.  At  the  very  beginning,  in  the  region  of  0(a), 
the  flow  develops  as  described  above.  However,  the  centrifugal  force 
becomes  one  of  the  dominant  forces  as  the  fluid  moves  into  the  next 
region,  a distance  of  0(a//a).  Much  of  the  flow  development  occurs 
within  a distance  of  The  variation  of  the  centrifugal  force 

and  the  pressure  gradient  directed  away  from  the  center  of  curvature 
on  a cross  section  of  the  pipe  is  small,  since  a is  generally  smaller 
than  one  for  most  curved  pipes  of  interest . The  resultant  of  the  nearly 


-5- 


uniformly  distributed  pressure  gradient  and  the  strong  centrifugal 
force  will  accelerate  the  fluid  in  the  central  core.  The  flow  is 
nealy  uniform  from  the  inner  to  the  outer  wall  of  the  tube  (away 
from  the  center  of  curvature) . The  boundary  layer  acts  as  a reser- 
voir, receiving  the  fluid  moving  towards  the  outer  wall,  and  also  as 
a source  of  fluid  leaving  it  at  the  inner  wall.  The  resulting  cross 
flow  forms  a stagnation-like  flow  locally  along  the  outer  wall  of 
the  pipe.  The  convective  effect  of  this  locally  stagnant  flow  pre- 
vents the  secondary  boundary  layer  from  growing.  Thus,  the  secondary 
boundary  layer  will  remain  thin  as  the  flow  asymptotically  approaches 
the  fully  developed  state  (Baura's  flow)  in  the  region  of  distance 
0(aRe/yl5).  There  are  three  length  scales:  a,  a/i/a,  aRe//D  for  this 

case.  The  solutions  of  regions,  distances  0(a/i/5)  and  0(aRe//D)  were 
given  by  Yao  and  Berger  (1975).  It  is  interesting  to  note  that  the 
size  of  the  region  where  Singh's  flow  is  valid  depends  on  the  values 
of  a.  It  shrinks  as  the  values  of  a increase. 

3.  When  both  a and  D are  large,  the  region  governed  by  Singh's 
flow  is  shrunk  to  zero.  (It  should  be  noted  that  when  we  say  that  a 
is  large  it  does  not  imply  that  the  value  of  a is  larger  than  one. 

It  simply  means  that  the  value  of  a is  not  close  to  zero.)  For  this 
case,  the  developing  flow  is  governed  by  the  solution  presented  by 
Yao  and  Berger  (1975) . 

HEATED  STRAIGHT  PIPE 

The  developing  flow  in  heated  straight  pipes  can  be  categorized 
into  cases  similar  to  the  developing  flow  in  curved  pipes.  It  should 
be  noted  that  the  fully  developed  flow  in  heated  straight  pipes  under 
the  condition  of  constant  wall  temperature  is  simply  the  Poiseuille 
flow,  and  the  ones  under  the  condition  of  constant  wall  heat  flux  are 
given  by  Morton  (1958)  if  ReRa  is  small  and  by  Mori  and  Futagami  (1967) 
if  ReRa  is  large.  Three  categories  of  the  developing  flows  are: 

(i)  When  e and  Gr  are  small,  the  entry  flow  deviates  slightly 
from  the  developing  Poiseuille  flow.  There  are  two 
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regions : 0(a)  and  O(aRe),  which  are  similar  to  the  entry 

flow  in  an  unheated  pipe. 

(ii)  When  e is  small  and  Gr  is  large,  the  axial-length  scales 

are  0(a)  and  0 (aRe^^^/Gr^^^)  for  the  constant  wall  heat 

1/2 

flux  and  0(a),  0(aRe/Gr  ) and  O(aRe)  for  constant  wall 
temperature.  The  secondary  flow  becomes  one  of  the  domi- 
nant flow  components  as  the  fluid  moves  into  the  region  of 
1/2 

0(aRe/Gr  ) for  constant  wall  temperature).  For  constant 
heat  flux,  the  flow  approaches  its  fully  developed  state 
at  the  end  of  region  0 (aRe^^^/Gr^^^) . For  the  case  of 
constant  wall  temperature,  the  temperature  difference  be- 
tween the  wall  and  the  core  flow  decreases  in  the  region 
of  O(aRe),  and  the  flow  approaches  its  fully-developed 
state,  Poiseuille  flow. 

(iii)  When  e is  large,  the  region  of  0(a)  does  not  exist.  The 
rest  of  the  flow  development  is  similar  to  case  (ii) . 

The  solution  presented  in  this  report  is  limited  to  the  region  distance 

■k 

0(a)  from  the  inlet  for  both  constant  wall  temperature  and  wall  heat 
flux  conditions  The  solution  exists,  apparently,  when  e is  small. 

As  the  fluid  enters  the  pipe,  a boundary  layer  develops  on  the 
wall  with  a secondary  flow  developed  inside  the  boundary  layer  by  the 
buoyancy  forces  due  to  the  temperature  difference  between  the  wall  and 
the  inlet  fluid.  The  central  inviscld  flow  will  be  accelerated  by  the 
displacement  of  the  axial  boundary  layer  like  that  in  an  unheated 
straight  pipe.  This  motion  generates  node-like  sinks  along  the  center 
line  of  the  pipe  and  causes  fluid  particles  to  move  toward  the  center 
of  the  pipe.  Simultaneously,  the  displacement  of  the  secondary  bound- 
ary layer  induces  a uniform,  downward  stream  on  the  cross-section  of 
the  pipe.  The  combination  of  the  downward  stream  with  the  radial 
stream  responds  to  the  development  of  two  vortices  in  the  central  core 
of  the  pipe.  A saddle-polnt-like  stagnation  point,  initially  coinci- 
dent with  the  center  line,  moves  downward  along  the  vertical  symmetry 


The  solution  in  the  region  of  0(a)  has  been  proved  to  be  valid 
when  z ^ 0.46/t^  for  external  flow. 
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llne  of  the  pipe,  and  the  saddle-point-like  stagnation  point  will 
stay  in  the  core  of  the  pipe  within  the  region  of  0(a).  The  axial 
velocity  profile  turns  counterclockwise  to  accommodate  the  downward 
stream  generated  by  the  displacement  of  the  secondary  boundary  layer. 

A favorable  pressure  gradient  develops  along  the  bottom  wall  of  the 
pipe;  on  the  contrary,  an  unfavorable  pressure  gradient  develops  along 
the  top  wall  of  the  pipe.  It  is  found  that  the  temperature  in  the 
central  core  of  the  pipe  does  not  change  within  the  region  of  0(a)  and 
is  identical  to  the  inlet  temperature.  It  is  expected  that  the  strength 
of  the  sinks  along  the  pipe  center  line  will  gradually  decrease  and  the 
downward-stream  pattern  will  become  dominant  as  the  fluid  moves  down- 
stream of  the  region,  a distance  0(a)  from  the  inlet  when  Gr  is  large. 

The  theory  developed  above  has  been  applied  to  the  "Colorado 
tube  test."  The  range  that  heating  stabilizes  the  boundary  layer  can 
be  found  in  the  region,  a distance  0(a)  from  the  inlet,  and  the  second- 
ary flow  will  be  negligible  in  this  region  since  the  values  of  £ for 
the  tests  are  very  small.  Downstream  of  the  region  of  0(a),  the  un- 
favorable pressure  gradient  along  the  top  wall  of  the  pipe  can  force 
the  flow  to  separate.  The  test  would  be  successful  if  the  flow  in  the 
test  section  can  be  restricted  to  the  region  distance  0(a)  from  the 
inlet.  The  size  of  this  region  shrinks  when  the  value  of  e increases. 

A detailed  evaluation  of  the  possible  test  range  of  the  laminar  flow 
is  given  in  Section  V. 
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II.  GOVERNING  EQUATIONS 


Near  the  inlet  it  is  natural  to  refer  lengths  to  the  radius  a of  the 

pipe,  and  velocities  to  the  uniform  inlet  velocity,  (Fig.  1).  The 

pressure  is  nondimensionallzed  by  p.  W.  where  p.  is  the  density.  The 

in  in  ^in 

dimensionless  temperature  is  defined  in  terms  of  the  inlet  temperature, 

T.  , and  the  constant  wall  temperature,  T,as0=(T-T.  )/(T  -T.  ), 

in  w in  w in 

or  the  constant  wall  heat  flux,  q , as  0 = (T  - T . ) k/R^/a  • q . The 

w in  w 

dimensionless  equations  of  motion  and  energy  with  the  Boussinesq  approx- 
imation in  cylindrical  coordinates  become 


1 3(rU)  , 1 W W ^ „ 
r 3r  r 34>  9z 


u|U  + V 3U 

3r  r 9(|)  9z  r 


9P 

-r h e (0  - 0 ) COS 

dr  c 


+ (v2u  _ IL  _ ^ 

Re  ^ 2 2 9*'’ 

r r 


9V  V 9V  , , 9V  . UV  3P  , .n  n ^ • a 

U — + — — + W — + — = — + e (0  - 0 ) sin 

dr  r 9(j)  9z  r rdcp  c 


+ (v2v  _ 

Re  2 2 9d)^ 

r r 


,,  3W  , V 9W  , ,,  aw 


9z  Re 


90  V 90  ^ ^ _J_  2 

9r  r 9(f)  9z  Pr  Re 


(2a) 

(2b) 

(2c) 

(2d) 

(2e) 


where 


2 2 2 

v2=^  + il_+l_^+^ 

. 2 ^ r 9r  2 ^ 2 ^ 2 

9r  r 9(J)  9z 


(3) 
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GravitaHon 


Fig.  1 — Coordinate  systems 

is  the  Laplace  operator  and  0^  is  the  temperature  along  the  center 
line  of  the  pipe.  The  parameters  e.  Re  and  Pr  are  defined  in  Eq.  (1) . 

The  entry  condition  is  uniform  inlet  axial  velocity  and  tempera- 
ture; the  reference  pressure  at  the  inlet  is  set  equal  to  zero.  It 
may  he  noted  that,  in  the  absence  of  viscosity,  the  exact  solution 
of  Eqs . (2)  satisfying  the  above  inlet  condition  and  slip  wall  con- 

dition is 


W = 1, 


U = V = P = 0=  O 


(4) 
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III.  SOLUTIONS 


As  the  fluid  flows  into  the  pipe,  the  viscous  forces  are  confined 
to  the  thin  boundary  layer  near  the  wall  of  the  pipe.  For  a heated 
pipe  the  secondary  flow  is  created  by  the  buoyancy  forces  inside  the 
thermal  boundary  layer.  The  ratio  of  the  thicknesses  of  the  thermal 
boundary  layer  to  the  momentum  boundary  layer  depends  on  the  Prandtl 
number.  Viscous  forces  and  heat  conduction  away  from  the  boundary 
layer  can  be  ignored;  the  flow  is  isothermal  and  inviscid  in  the  central 
core  of  the  pipe.  The  core  flow  is  accelerated  due  to  the  displacement 
effect  of  the  boundary  layer  and  fluid  particles  will  be  pushed  from 
the  wall  toward  the  center  of  the  pipe.  Simultaneously,  a downward 
stream  is  developed  due  to  the  displacement  effect  of  the  secondary 
boundary  layer.  The  combination  of  the  radial— direction  motion  with 
the  downward  one  gives  the  stream  pattern  of  two  developing  vortices 
on  the  cross-section  normal  to  the  axis  of  the  pipe.  The  downward 
stream  will  also  skew  the  axial  velocity  by  accelerating  it  below  the 
center  and  decelerating  it  above  the  center  of  the  pipe.  The  analysis 
shows  that  the  development  of  the  secondary  flow  due  to  the  heating 
effect  near  the  entrance  when  e is  small  can  be  obtained  by  perturbing 
the  solution  of  the  developing  flow  in  an  unheated  pipe. 

ZEROTH-ORDER  SOLUTION  IN  THE  INVISCID  CORE 

The  solution  of  the  inviscid  core  flow  can  be  obtained  by  expanding 
the  dependent  variables  into  series  in  6 which  are 


u = u 

+ 

SU^ 

4-  ... 

o 

> 

II 

> 

+ 

5V^ 

+ • . • 

w = w 

+ 

6W^ 

+ ... 

P = Po 

+ 

6Pi 

+ • . . 

0 = 00 

+ 

60^ 

+ . . • 

>c  = ®S 

+ 

60^ 

+ • . . 

(5) 
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where  6,  the  expansion  parameter,  can  be  determined  by  matching  with  the 
zeroth  order  boundary-layer  flow.  It  can  be  shown  that  it  is  equal  to 
1//Re.  The  governing  equations  of  the  first  order  can  be  obtained  by  tak- 
ing the  limit  of  Re  °°  from  Eqs . (2).  The  zeroth  solutions  which  satisfy 
the  uniform  inlet  velocity  and  temperature  conditions  and  the  slip  condi- 
tion at  the  pipe  wall  are  the  undisturbed  flow,  which  is  given  in  Eq.  (4). 

ZEROTH-ORDER  BOUNDARY-LAYER  FLOW 

Near  the  pipe  wall,  the  viscous  forces  and  heat  conduction  normal  to 
the  wall  become  important.  The  radial  coordinate  r is  stretched  to  reflect 
this  physical  fact.  Accordingly,  we  introduce  the  inner  variables,  as  in 
the  classical  boundary  layer, 

r = 1 - 6r,  U = -5u,  V(r,c|),z)  = v(r,()),z),  etc.  (6) 

A combination  of  Eqs.  (2)  and  (6),  after  neglecting  the  smaller  order  terms, 
yields 


9u  3;v  ^ _ 

97  ^ 


2 

9v  , 9v  , 9v  9P  , Q \ • A 3 

- 94.  9z  94.  c g-2 


9w  , 9w  , 9w  9P  , 3 w 

+ + ^ 97  = “■^  + 7=2 

9r  9r 


(7) 


99  99  , 99  1 9 9 

u + v — + w — = :r r 

97  972 


The  analysis  of  the  first-order  boundary-layer  flow  can  be  easily  extended 
to  the  case  in  which  fluid  properties  depend  on  the  temperature.  For  most 
liquids,  the  principal  departure  from  constant-property  flow  is  due  to 
the  viscosity  variation.  Including  the  variable  viscosity  of  the  fluid, 
the  terms  of  the  viscous  forces  in  Eqs.  (7)  will  be  replaced  by 
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8 

->  — 


dr^ 


3^w 

37^ 


(-S) 

ir  \ 3r  / 


3r 

3 


(8) 


where 


N = Ji- 


is  the  viscosity  normalized  by  the  viscosity  in  the  central  core.  The 
associated  boundary  conditions  of  Eqs.  (7)  are 


Atz  = 0,u=v=e  = 0,w  = l 
Atr=0,  u=v=w=0, 
e = 1,  or 


0-  = -1 
r 


Asr-^°°,v->0,w->l,  0->O 


(entrance  condition)  ^ 
(no-slip  condition  at  wall)  j 
(constant  wall  temperature).^ 
(constant  wall  heat  flux) 

(matching  condition  with 
zeroth-order  inviscid 


(9) 


Under  the  constant  wall  temperature,  the  solutions  of  Eqs.  (7) 
satisfying  the  conditions  (8)  have  been  published  by  Yao  and  Catton 
(1976a, b)  for  Pr  = 0.01,  1,  10  of  the  constant  property  flow  and  Pr  = 8 
of  the  water  flow  (variable  viscosity).  Their  results  will  be  sum- 
marized briefly  below  in  order  to  make  this  report  self-contained. 

The  dependent  variables  can  be  expressed  as 


w = fp'  (p)  + e(2z)  F^'  (p)  • cos  (t>  + . 


V = e(2z)  • F2' (P)  ' sin  <j)  + . . . 

u = -^  (pf„'  - f„)  + e(2z)^/^pF  ' - 5F  - F ) 
0 0 112 


cos  tb  + 


0 = s(2z)  • G(p)  • cos  <f)  + 


(10a) 

(10b) 

(10c) 


(lOd) 
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where  n = r//^  is  the  Blasius  variable.  Substitution  of  Eqs.  (10) 
into  Eqs.  (7)  and  collecting  terms  of  equal  order  of  £ yields 


(N  f 
^00 


+ ^o^o"  = 0 


+ Pr  foG^'  = 0 


(lla) 

(llb) 


and 


(NqF^")’  + foF^"  - 4fo'F^’  - 5fo”F^  + V’F^  = (A  • G • fQ-) 


W - 2fo’^2'  = -®0 


(12) 


where 


|^G"+fQG'  - 4fo’G=  -eo'(5F^+F2) 


N 


0 ■ 1 + «(T„  - 


A = 


- bn> 

[i  + - b„)»o]‘ 


(13) 


and  a is  the  coefficient  of  the  variable  viscosity.  For  water  in  the 
temperature  range  between  40°F  (4.4°C)  and  100  F (37.8  C) , a ^ 0.0151  ( F) 
for  the  constant  property  flow,  a = 0.  The  boundary  conditions  (9) 
become 


^o' 

= 0 

and 

1 

at 

q = 0 

^o' 

-5-  0 

and 

o 

CD 

0 

as 

P ->■  00 

(14a) 
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and 


= F^'  = G = 0 at  n = 0 | 

Fi',  F2',  and  G->'^  as  n-^“  i 


(14b) 


The  numerical  values  of  stream  functions  f^,  0^,  F^,  F^,  and  G can  be 
found  in  Yao  and  Catton  (1976a, b). 

The  velocity  normal  to  the  wall  along  the  outer  edge  of  the  bound- 
ary layer  is  required  to  solve  the  second-order  inviscid  core  flow 
(matching  conditions) . The  matching  conditions  can  be  obtained  by 
taking  the  limit  of  Eq.  (10c)  for  n which  gives 

U (r  = 1)  = - limit  5 -i-  (qf  ' - f ) + 

X (nF^*  - 5F^'  - F^)  cos  (])  + ...| 


P-j  'X  / 0 

e3t(2z)  cos  (()  + — (15a) 


where 


3 = limit  (n  - f„)  (15b) 

1 n ^ 

g = limit  (5F  + F ) (15c) 

^ n ^ ^ 


The  values  of  and  62  =i0.01,  1.0,  10  of  the  constant  property 

flow  and  for  Pr  = 8 of  the  water  flow  are  taken  from  Yao  and  Catton 
(1976a, b)  and  are  listed  in  Table  1. 
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Table  1 

CONSTANTS  FROM  THE  COMPUTATION  OF  THE  BOUNDARY -LAYER  FLOW 


A.  Constant  T 

w 


Constant 

Property^ 

Variable  Viscosity  (Pr 

= 8)t> 

Pr 

0.01 

1.0 

10.0 

n^at 

0 

0.5 

1 

^1 

1.2167 

1.2167 

1.2167 

^1 

1.2167 

1.0617 

0.9497 

®2 

3.9454 

0.8534 

0.3158 

0.3522 

0.3106 

0.2835 

B.  Constant  q (constant  property) 


Pr 

bX 

0.01 

1.0 

10.0 

^1 

1.2167 

1.2167 

1.2167 

B3 

23.3511 

0.9062 

0.1524 

i 

^Yao  and  Catton,  1976a. 
^Yao  and  Catton,  1976b. 


For  the  case  of  constant  wall  heat  flux,  the  dependent  variables 
can  be  expanded  as 


w = fg'  + • F^'  • cos  (f)  + ...  (16a) 

V = e(2z)^^^  • F^'  • sin  c(>  + ...  (16b) 

u = (nf„'  - f„)  + e(2z)^  • (nF  ' - 6F  - F ) • cos  cf)  + . . . (16c) 

0 0 11/ 

0 = -fli.  • 0Q  + e(2z)^  • 


G • cos  (j)  + . . . 


(16d) 
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The  governing  equations  for  £q,  6q,  F^,  F^,  and  G are 


f "'  + f f " = 0 
0 0 0 


^o"  + Vo’  - ^o’^o  = ° 


F + f F " - 5f  'F  ' + 6f  "F  = -f  "F 
1 0 1 0 1 0 1 0 2 


F "'  + f F " - 3f  'F  ' = -9 
2 0 2 0 2 0 


^G"+foG-  - 6fQ’G=-9QF^'  - 9q'(6F^_F2) 


(17a) 

(17b) 

(17c) 

(17d) 

(17e) 


The  boundary  conditions  of  Eqs.  (18)  for  constant  wall  heat  flux  are 


‘o  - - 0-  ®o’  - “ 


^1,  ®o  ^ ^ 


(18a) 


and 


Fi  = Fi'  = F2  = F2'  = G'  = 0 at 
F^',  F^',  G ->  0 as 

The  normal  velocity  along  the  outer  edge  of  the  boundary  layer 
which  is  the  matching  condition  for  the  second-order  invlscld  core  flow 
is 


n = 0 


(18b) 


Pi  2 

1-  * £ (2z)  • cos  (f)  + . . . 


U^(r  = 1)  = - 


(19a) 
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where 


3 


3 


limit 

n ->  oo 


(6F^  + F^) 


is  also  listed  in  Table  1. 
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FIRST-ORDER  INVISCID  CORE  FLOW 

The  equations  of  the  first-order  inviscid  core  flow  are  found  by 
substituting  Eqs.  (4)  and  (5)  into  Eqs.  (2).  They  are 


1 

„ 1 + r. 

= 0 

(20a) 

r 

3r  r 3(})  3z 

au- 

1 _ 

- 3r  ^(®1  - ®1  > 

cos  <f) 

(20b) 

3P 

1 

8z 

r 3(J)  11 

sin  (j) 

(20c) 

8W, 

3P- 

1 

1 

(20d) 

3z 

3z 

30, 

1 _ 

0 

(20e) 

3z 

The  temperature  of  the  inviscid  core  flow,  up  to  this  order,  is 
still  not  changed  and  is  equal  to  the  inlet  temperature.  This  is  con- 
firmed by  Eq.  (20e) , the  solution  of  which  is  0^  = 0.  Therefore,  there 
is  no  buoyancy  force  acting  on  the  fluid  in  the  central  core  of  the  pipe. 
Equations  (20a)  through  (20d)  can  be  separated  into  two  parts:  the  accel- 

erating axial  flow  due  to  the  displacement  effect  of  axial  boundary  layer 
and  the  downward  stream  due  to  the  displacement  effect  of  the  secondary 
boundary  layer . The  equations  governing  the  two  phenomena  are  found  by 
substituting 
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^1 

= ^10-^ 

^''ll 

COS 

+ • . . 

^1 

= EV^^  , 

sin  (j) 

+ . . 

= ^10^ 

^^11 

cos 

4> 

+ . . . 

eP^I 

cos 

<f> 

+ • . . 

into  Eqs.  (20).  After  collecting  the  terms  without  e and  with  e,  we 

, 9(rU,,)  3W, 


r 3r 


^10^  + i;io  = 0 


8z 


^^10  ^^10 

9z  9r 


=>“io  ®10 


3z 


3z 


and 


1 3 V ^ ^^11 

i = 0 

r 3r  r 3z 


3U 


11 


3P 


11 


3z 


3r 


'\l  "ll 


3z 


3W, 


11 


3P 


11 


3z 


3z 


Displacement  Effect  (o(5)) 

The  appropriate  matching  condition  for  Eqs.  (22)  is 


(21) 


have 

(22a) 


(22b) 


(22c) 


(23a) 


(23b) 


(23c) 


(23d) 
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Uio  = - 6j^/(2z)  ^ at  r = 1 


(24a) 


and  the  entry  conditions  are 


ho  " “lO  - “ 


at  z = 0 


(24b) 


Integrating  Eq.  (22c)  with  respect  to  z and  using  (24a)  gives 


W = - P 
10  ^10 


(25) 


Eliminating  and  from  Eqs.  (22a),  (22b)  and  (25)  results  in 


’"ho  . 1 ’ho  . ’"ho  _ „ 


8r 


r 9r 


(26) 


3z 


The  solution  of  Eq.  (26)  satisfying  condition  (24a)  can  be  found  by 
Fourier  Transform  in  the  sense  of  generalized  functions  (see  Lighthill, 
1958).  It  is 


10 


-/iT  Jo  “ 


Io^(«r) 


(a) 


sin  az  da 


where  I's  are  modified  Bessel  functions. 

Substituting  Eq.  (27a)  into  Eqs.  (22b)  and  (22c),  we  obtain 


W 


10 


1 

TT 


a ^ . V — sin  az  da 

Ij^(a) 


(27a) 


(27b) 


and 

/*“  h Ii  (ar) 

U i I a ^ (cos  az  - 1)  da 

“ /IT  70  h<“> 


(27c) 
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The  asymptotic  values  of  Eqs.  (27)  for  large  z can  be  easily  found: 


""lo  = - ^0  = ®1 


J-  r 

2(2z)^  • 1 + Js  2 

L (2z)^ 


+ • • c 


(28a) 


• r • (2z) 


(28b) 


Equations  (27)  or  (28)  represent  the  accelerating  flow  due  to  the  displace 
ment  effect  of  the  boundary  layer  in  an  unheated  pipe. 


Secondary  Displacement  Effect  0(5e) 

For  the  constant  wall  temperature,  the  matching  condition  for  Eqs. 

(23)  is  found  from  Eqs.  (15)  and  is 

Uii  = B2(2z)^^^  at  r = 1 (29a) 


The  entry  conditions  are 


P = W,  = 0 at  z = 0 
11  11 


(29b) 


Solving  Eqs.  (23)  with  conditions  (29)  yields 


6B. 


/F  •'o 


r“  -5/2 

Jn  “ I 

•'n  I 


I^(ar) 


^(a)  + 12(a) 


sin  az  • da 


(30a) 


68. 


^11  Jo  “ ^0^“^  ^ ^2^“^ 


, I,  (ar) 

5/2  1 — (cos  az  - 1)  da  (30b) 


and 


11 


369  s/9  3-  l2(ar) 

2 ‘ ° (cos  az  - 1)  da 


I 


Iq(«)  + 12(a) 


(30c) 
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For  a large  z , they  can  be  approximated  by 


W 


11 


11 


3B. 


2 

r 

4 


(2z)-3/2  + 


(31a) 


and 


-Vii  = ~ B2(2z)^^^  + . . . (31b) 

For  the  case  of  constant  wall  heat  flux,  the  matching  condition  is, 
from  Eqs.  (19) , 


Uii  = 6^(22)  at  r = 1 


With  conditions  (29b)  and  (32),  the  solutions  of  Eqs.  (23)  are 


I^(ar) 

„ “■  • I„(a)  + I^(a)  • 


46^  r(2z) 


Vfi  = - 63(2z)^ 


Uii  = 63(2^)' 


(32) 


(33a) 

(33b) 

(33c) 


where  5^ (a)  is  the  first  derivative  of  the  delta  function. 

The  analysis  can  be  extended  systematically  to  obtain  higher  order 
terms  of  the  solution  in  the  boundary  layer  and  the  core.  The  temperature 
of  the  center  core  is  not  changed  up  to  0(o£);  however,  the  downward  convec- 
tion, 0(6e),  will  carry  the  hotter  fluid  from  the  boundary  (-  ir/2  ^ cj)  ^ it/2) 
into  the  cooler  core  and  will  increase  the  temperature  of  the  core  flow. 

The  small  temperature  variation  is  not  considered  in  this  report. 
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IV.  RESULTS  AND  DISCUSSION 


BOUNDARY-LAYER  FLOW 

A detailed  discussion  of  the  boundary-layer  flow  for  the  case  of 
constant  wall  temperature  can  be  found  in  Yao  and  Catton  (1976a, b)  and 
will  not  be  repeated  here. 

For  the  case  of  constant  wall  heat  flux,  the  numerical  results 

of  functions  0^ , , F^,  F^',  F^  and  F2'  are  plotted  in  Figs.  2-7.  The 

values  of  0 for  the  case  of  constant  q are  larger  than  those  for  the 
0 "w 

case  of  constant  T . The  wall  temperature  for  the  case  of  constant  q 
w 

increases  downstream,  and  its  values  can  be  calculated  from  Eq . (16d) . 

The  values  of  0^  and  for  Pr  = 0.01,  1.0,  and  10  are  listed  in  Table 

2.  In  Fig.  3 the  values  of  the  perturbed  temperature  are  positive  near 

the  wall  and  negative  near  the  core  for  the  case  of  constant  q^.  For 

the  case  of  constant  T they  are  always  negative  and  are  not  plotted, 

w 

their  values  are  too  small  to  be  shown  in  the  scale  of  Fig.  3. 

The  stream-function  F^  shown  in  Fig.  4 indicates  that  their  values 
are  much  larger  than  that  for  the  case  of  constant  T^.  The  induced 
axial  velocity  profile  has  more  profound  effects  for  the  case  of  con- 
stant q . Similarly,  the  buoyancy  force  induces  stronger  secondary 
w 

flow  at  the  constant  wall  heat  flux  than  at  the  constant  wall  temper- 
ature shown  in  Figs . 6 and  7 . 


Table  2 

COEFFICIENTS  FOR  WALL  TEMPERATURE  AND  SHEAR  STRESS 


Pr 

0o(O) 

G(0) 

F^"(0;/fQ"(0) 

0.01 

9.1173 

4.6807 

1.2085 

7.1801 

1.0 

1.5409 

1.0053 

0.0990 

0.7914 

10.0 

0.7087 

0.0958 

0.0198 

0.2309 
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Fig.6 — Function  F2 
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SHEAR  STRESS 

The  local  shear  stress  at  the  wall  can  be  computed  from  the  equation 


T 

rz 


r=l 


and 


r(J) 


3v 

9r 


r=l 


Introducing  the  series  expansions  (16) , the  relative  importance  of  the 
secondary  flow  on  the  axial  shear  stress  can  be  found  from 


T 

rz 


<Rz>o 


1 + e(2z) 


5/2 


F^"(0) 

VW 


cos  4>  + . . . 


(34) 


The  circumferential  shear  stress  can  be  shown  to  be  proportional  to 


r(j) 


e(2z) 


3/2 


F2"(0)  • sin  (j) 


(35) 


Values  of  F^"  (0) /f q" (0)  and  F^'HO)  are  given  in  Table  2.  Equations  (16d) 

and  (34)  indicate  that  the  secondary  flow  effect  on  the  wall  temperature 

5/2 

and  shear  stress  grows  rapidly  downstream,  and  is  proportional  to  z 
This  means  that  an  initially  small  secondary  effect,  which  may  be  treated 
as  a second-order  effect  in  the  region  a distance  0(a)  from  the  inlet, 
becomes  a dominant  flow  component  farther  downstream. 

INVISCID  CORE  FLOW 

Equations  (27a, b)  and  (30a)  are  numerically  integrated  by  using  the 
Fast  Fourier  Transform.  The  details  of  the  numerical  method  are  described 

in  Appendix  A.  "^lO'^^l^  '^11'^ ^2  ~^ll'^^2'^  evaluated 

at  r = 0.5  as  functions  of  z,  and  plotted  in  Fig.  8.  The  values  of  W^g/3^ 
(or  at  r = 0 differ  only  slightly  from  their  values  at  r = 0.5 

and  are  not  present  in  the  figure.  Along  the  center  line,  functions  and 

are  identical  to  zero.  is  the  flow  acceleration  due  to  the  dis- 

placement effect  of  the  axial  boundary  layer.  The  flow  is  quickly  acceler- 
ated near  the  entrance  and  the  acceleration  declines  and  is  proportional  to 
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Fig.  8 — Displacement  effect 


Pi 


or 


Pi 


and 


secondary  flow  effect 


p2 


or 


iM 

p2 
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downstream.  The  distributions  of  (or  P^q)  on  the  cross-section 

normal  to  the  pipe  axis  are  given  in  Fig.  9.  Near  the  entrance  z = 0.1, 
the  axial  velocity  is  higher  close  to  the  boundary  layer  than  along  the 
center  line  of  the  pipe.  The  delay  of  the  flow  acceleration  near  the 
center  line  causes  the  axial  velocity  profile  to  become  concave.  The 
velocity  profile  has  its  maximum  velocity  off  the  center  line.  The  peaks 
are  eroded  downstream  and  eventually  disappear  at  z = 2.5;  see  Fig.  9. 

A similar  interpretation  can  be  given  to  k-|^Q5  which  is  the  pressure  drop. 

The  displacement  effects  of  the  secondary  boundary  layer  are  repre- 
sented by  functions  ^11'  Physically,  the  secondary  boundary 

layer  causes  the  axial  flow  to  turn  counterclockwise  around  the  horizontal 
line  passing  through  the  center  of  the  pipe  and  normal  to  the  axis  of  the 
pipe,  Y = r cos  ())  = 0,  and  does  not  cause  any  net  mean  acceleration  of  the 
flow.  The  turning  rate  is  higher  near  the  entrance  and  gradually  matches 

-h 

the  rate  which  is  proportional  to  z 

For  the  case  of  content  T , the  asymptotic  expansions,  Eqs . (28a)  and 

w 

(31a),  of  Eqs.  (27a, b)  and  (30a)  for  large  z's  are  also  plotted  in  Fig.  8. 
They  show  that  the  asymptotic  values  of  (or  P^^q)  match  their  exact 

value  approximately  at  z = 2.  For  (or  asymptotic  values 

start  to  be  valid  at  z = 4.  The  comparison  of  the  asymptotic  expansions 
with  the  exact  evaluation  of  the  sine  Integrals  suggests  that  simpler 
asymptotic  expansions  can  be  used  for  practical  application  when  z is 
larger  than  four,  which  will  cover  most  ranges  of  practical  Interest  when 
e is  small.  The  velocity  components  and  the  pressure,  in  terms  of  their 
asymptotic  forms,  are 


P 


+ e • 362  • ^ cos  (f)  • 


w = 1 - 

U = -5 
V = -5e6 


+ . 


sin  (j) 


? -2 
r (2z) 


[1  + ha  - ^)(2z)  ^ + 


1 r 3/2 

. . - 062  (2z)  • cos 


c|)  + . 


(36a) 

(36b) 

(36c) 


(36d) 
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For  the  case  of  constant  q , the  terms  of  0(6e)  in  Eqs.  (36)  are  replaced 

w 

by  Eqs . (33) . 

PRESSURE  DROP 

The  pressure  head  will  drop  along  the  pipe  as  the  flow  is  accelerated. 
The  first  term  of  Eq.  (36a)  represents  the  pressure  drop  in  an  unheated 
pipe  due  to  the  displacement  of  the  axial  boundary  layer.  The  pressure 
drops  faster  close  to  the  edge  of  the  boundary  layer  than  along  the  center 
line  of  the  pipe  where  r = 0.  The  difference  in  the  pressure  distribution 
on  the  cross-section  of  the  pipe  disappears  for  large  z,  and  the  pressure 
distribution  approaches 


P ~ -6  • 2B^(2z)^ 


(37) 


The  pressure  distribution  is  disturbed  when  the  pipe  wall  is  heated,  which 

is  represented  by  the  second  term  of  Eq.  (36a).  This  term  can  be  rewritten 
^ 

in  (X,Y,z)  coordinates  as  6e  3$2Y(2z)^  for  constant  T^,  and  as  6e  86^  Yz 
for  constant  q^,  which  shows  that  the  displacement  effect  of  the  secondary 
boundary  layer  introduces  an  unfavorable  pressure  gradient  over  the  upper 
half  of  the  pipe  flow,  Y > 0,  and  a favorable  one  over  the  lower  half  of 
the  pipe,  Y < 0. 


AXIAL-VELOCITY  PROFILES  OF  THE  CORE  FLOW 

The  axial  velocity  of  the  core  flow  can  be  viewed  as  a superposition  of 
these  components.  The  first  one  is  the  undisturbed  flow  at  the  inlet,  i.e., 
W = 1.  The  second  one  is  the  accelerated  flow  due  to  the  displacement 
effect  of  the  axial  boundary  layer  in  an  unheated  pipe.  The  last  com- 
ponent is  due  to  the  displacement  of  the  secondary  boundary  layer.  Similar 
to  the  pressure  distribution,  the  third  term  in  (X,Y,z)  coordinates  is 

-6e36.,Y(2z)  ^ for  constant  T , and  -5e  86  Yz  for  constant  q , which  means 
2 w w 

that  the  axial  velocity  profile  turns  counterclockwise  around  Y = 0.  In 

the  region  of  z ~ 0(a),  both  the  second  and  the  third  terms  increase  as 

for  constant  T , and  as  z for  constant  q . However,  the  second  term 
w w 

may  be  dropped  when  the  third  term  is  still  growing  downstream  beyond 
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the  region  of  0(a)  when  Gr  is  not  small.  For  the  cases  of  small  Gr  (or 
ReGr) , the  second  and  the  third  terms  will  continue  growing  and  the  flow 
will  approach  its  fully  developed  state. 

VELOCITIES  ON  THE  CROSS-SECTION 

The  radial  velocity  component  U is  given  in  Eq.  (36c).  The  first 

term  is  due  to  the  displacement  effect  of  the  axial  boundary  layer,  which 

displaces  the  fluid  away  from  the  wall  and  toward  the  center  of  the  pipe. 

The  second  term  represents  the  motion  when  the  fluid  leaves  the  boundary 

layer  over  the  upper  half  of  the  pipe  (tt/2  < (j)  < 3tt/2),  and  enters  the 

boundary  layer  over  the  lower  half  of  the  pipe  (-  tt/2  < (j)  < tt/2).  The  first 

kind  of  motion  declines  downstream,  and  the  second  one  increases  due  to 

heating.  They  become  the  same  order  of  magnitude  when  z ~ 0(a/y^)  (or 

0(a/e  ^ )).  Beyond  this  point,  the  analysis  is  no  longer  valid.  For  the 

case  of  large  Gr,  the  second  term  becomes  dominant  and  the  axial  length 

2/5 

scale  will  be  a/v^  (or  a/e  ' for  the  constant  q case).  For  a small  Gr , 

w 

2/5 

a/v^  (or  a/e  ' for  the  constant  q^  case)  is  larger  than  a • Re.  This 
suggests  that  the  next  length  scale  will  be  a • Re.  Physically,  it  suggests 
that  the  flow  development  will  be  similar  to  the  unheated,  straight  pipe  and 
the  solution  of  the  slightly  heated  pipe  can  be  obtained  by  perturbing  the 
solution  of  an  unheated  pipe.  Along  the  line  of  (j)  = 0,  where  V = 0,  there 
is  a stagnation  point  at 


r 


(contant  T ) 
w 


(constant  q ) 
w 


(38) 


This  is  also  a saddle  point,  point  P in  Fig.  10,  which  is  similar  to  the 
entry  flow  of  a curved  pipe  found  by  Singh  (1974).  In  a heated  pipe, 

(or  eB^/Bj^)  is  smaller  than  one.  Table  1,  and  the  analysis  requires 
the  quantity  [e(2z)^]  to  be  smaller  than  one.  Therefore,  the  saddle  point  P 
can  never  move  into  the  boundary  layer  within  the  region  of  0(a). 
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Equation  (36d)  shows  that  the  circumferential  velocity  has  a maximum 
at  c()  = tt/2  and  is  equal  to  zero  at  both  the  top  and  the  bottom  of  the 
pipe,  <p  = 0,TT. 


Top 


Pig  10 — Streamlines  on  the  cross  section 

STREAMLINES  ON  THE  CROSS-SECTION 

For  constant  T the  projection  of  the  streamline  on  a cross-section 
w 

can  be  computed  by 


3 r(2z)  - e3^  cos  (p(2z)^''^ 

1 dr  ^ ^ 2 (39) 

""  £32  sin  cj)(2z)^/^ 

Equation  (39)  shows  that  all  streamlines  are  tangent  to  the  vertical  line 
(})  = 0,TT,  and  are  plotted  in  Fig.  10.  If  z is  not  too  large,  the  second 
term  in  the  numerator  of  Eq.  (39)  can  be  neglected  and  Eq.  (39)  can  be 


integrated  to  give 


tan  Y = £ 


+ constant 


(40) 


Bit 


For  e = 0,  an  unheated  straight  pipe,  Eq.  (40)  becomes  (p  constant.  The 
streamlines  on  the  cross-section  are  a radial  straight  line.  Within  the 
limit  of  the  analysis,  the  second  term  on  the  numerator  of  Eq.  (39)  can  never 
be  larger  than  the  first  one.  Nevertheless,  stretching  our  interpretation 
of  Eq.  (39)  somewhat  and  considering  it  for  large  z can  predict  what  may 
occur  for  the  flow  development  farther  downstream.  Equation  (39)  becomes, 
after  neglecting  the  first  term  of  the  numerator  and  integrating, 

X = constant  (41) 


which  shows  that  the  streamlines  on  the  cross-section,  due  to  the  dis 
placement  effect  of  the  secondary  boundary  layer,  are  vertical  straight 
lines.  The  magnitude  of  this  downward  flow  is 

A/  + \l^  = 62(2z)^^^  (^2) 

which  increase  downstream.  This  flow,  when  Gr  is  large,  will  eventually 
become  the  dominant  flow  pattern,  which  is  observed  by  Mon  and  Futagami 
(1967).  Also,  the  downward  stream  forms  a stagnation-like  flow  locally 
along  the  bottom  wall  of  the  pipe.  The  convective  effect  of  this  locally 
stagnant  flow  prevents  the  boundary  layer  from  growing.  Thus,  the  boundary 
layer  will  remain  thin  as  the  flow  moves  downstream.  Therefore,  the 
flow  acceleration  due  to  the  displacement  effect  of  the  axial  boundary 
layer  will  fade  out.  This  suggests  that  the  flow  development  in  the  region 
of  0(a//e)  (or  0(a/e^''^))  is  mainly  in  changing  the  slope  of  the  axial 
velocity  profile  to  balance  the  gradually  enhanced  flow  in  the  secondary 
boundary  layer.  Also,  the  unfavorable  pressure  gradient  near  the  top  wall 
of  the  pipe  gradually  grows  and  eventually  triggers  the  separation  of  the 
boundary  layer  farther  downstream.  The  phenomenon  of  the  flow  development 
will  be  similar  to  that  of  a curved  pipe  given  by  Yao  and  Berger  (1975) . A 
similar  interpretation  can  be  made  for  the  case  of  constant  q^. 
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V.  THE  TUBE  TEST  IN  COLORADO  (Barker  and  Jennings,  1977) 


The  dimensions  of  the  tube  used  in  Colorado  are  shown  in  Fig.  11. 
The  heating  condition  is  assumed  to  be  constant  wall  temperature.  The 
range  of  Re,  Gr , and  e are  computed  on  the  basis  of  the  preliminary 
test  arrangement: 


T.  = 52°F 
in 

a = 0.0151  (“F)“^ 

Y = 0.49  X io~^  (°F)“^ 
a = 2"  = 0.167' 

g = 32.2  ft/sec^ 

V = 1.31  X 10“^  ft^/sec 
Pr  = 8 

W.  =5-30  ft/sec 
in 

AT=T  -T.  =0-  40° 

w in 

W.  a c £ 

Re  = = 1.27  X 10  - 1.02  X 10 

V 

Gr  = = 0 - 3.49  X 10® 


Gr  “4 

e = = 0 - 2.17  X 10  ^ 

Re^ 

6 = = 9.90  X io~^  - 2.81  X 10 


From  the  previous  section  we  know  that  the  secondary  flows  induced 
by  the  buoyancy  forces  are  0(6e)  in  the  core  flow  and  0(e)  in  the  boundary- 
layer  flow;  these  are  very  small  in  the  region  distance  0(a)  from  the 
entrance,  and  can  be  neglected.  The  test  data  obtained  in  the  region  can 
be  interpreted  as  being  equivalent  to  what  would  be  obtained  for  an  ex- 
terior flow,  except  that  there  is  a slight  acceleration  of  the  flow  due 
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4 


- 21  ft 


W jn  = 5 ~ 30  ft/ sec 

Tin=52°F 

Tw  ~^in  ~0~40°F 


Fig.  11  — Tube  test  in  Colorado 


to  the  displacement  of  the  axial  boundary  layer.  The  flow  acceleration 
can  be  calculated  from  Eq.  (36b),  which  is,  for  z = 4, 


AW  = 


W - W. 

xn 


W. 

xn 


6 


1^2(2z)^ 


(43) 


The  range  of  AW  for  21-ft-  and  45-ft-long  pipes  is  given  in  Table  3. 

Table  3 

FLOW  ACCELERATION 


(Pr  = 8) 

\^AW 

aAT 

21-ft  pipe  (%) 

45-ft  pipe  (%) 

0 

1.52  - 4.32 

2.29  - 6.49 

0.5 

1.33  - 3.77 

1.99  - 5.66 

1 

1.19  - 3.37 

1.78  - 5.61 

which  shows  that  the  flow  acceleration  decreases  when  the  heating  rate 

The  upper  bound  of  the  flow  acceleration  is  about  6.5  percent. 


increases . 
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In  general,  the  flow  acceleration  will  stabilize  the  boundary  layer;  the 
quantification  of  this  effect  is  being  studied  and  will  be  reported  when 
it  is  completed. 

In  view  of  the  above  examination,  the  data  obtained  from  the  tube 
test  will  be  valuable  in  simulating  external  flow  if  the  flow  region 
in  the  tube  falls  in  the  solution  of  the  region  distance  0(a)  from  the 
entrance.  However,  the  size  of  the  region  0(a),  the  solution  of  which 
is  presented  in  the  previous  sections,  is  determined  by  the  magnitude  of 
e.  From  a preliminary  investigation  of  the  next  region,  distance  0(aRe/Gr^) 
from  the  entrance,  it  is  indicated  that  the  size  of  the  region  0(a)  is 
about 


z 


0.46 

/e 


0.46W. 
xn 

/ygaAT 


W. 


28.36 


(44) 


This  shows  the  size  of  the  region  distance  0(a)  from  the  entrance  is 
about  7.5  ft  from  the  entrance  when  the  inlet  velocity  is  10  ft/sec 
and  At  = 40 °F  (condition  I) . This  distance  increases  when  the  inlet 
velocity  Increases  and  AT  decreases.  Its  value  equals  approximately 
44.9  ft  from  the  entrance  when  the  inlet  velocity  = 30  ft/sec  and 
AT  = 10°F  (condition  II).  Therefore,  under  condition  II,  the  flow  in 
the  test  tube  is  governed  by  the  solutions  presented  in  the  previous 
section  and  can  be  used  to  model  the  external  flow.  However,  under 
condition  I,  only  the  flow  of  the  first  7.5  ft  from  the  entrance  of 
the  test  tube  can  be  used  to  simulate  the  exterior  flow.  From  7.5  ft 
to  the  exit  of  the  tube,  the  flow  is  governed  by  the  solution  of  the 
region  distance  0(aRe/Gr^)  from  the  entrance,  since  the  values  of  Gr 
are  large  in  those  tests. 

From  the  tube  test  in  Colorado,  it  has  been  found  that  it  is  neces- 
sary to  increase  the  wall  temperature  and  the  flow  rate  simultaneously 
in  order  to  obtain  the  maximum  wall-heating  effect  on  the  flow  transi- 
tion. Otherwise,  the  measured  transition  Reynolds  number  was  inevitably 
much  lower.  This  is  apparent  because  the  "overheating"  causes  strong 
secondary  flow  and  destabilizes  the  boundary  layer.  This  phenomenon 
is  demonstrated  in  Fig.  12.  The  dot-circles  are  the  "successful"  data 
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from  the  tube  test,  and  the  two  straight  lines  indicate  the  possible 

location  of  the  end  of  the  region  of  0(a),  where  the  buoyancy  is  ex- 

2 

pected  to  be  negligible.  The  value  ex  = 0.21  is  the  line  beyond  which 
the  buoyancy  effect  becomes  important.  Physically,  the  slope  of  the 
line  can  be  viewed  as  the  "slope"  of  the  "laminar  path."  A preliminary 

r\ 

analysis  indicates  that  the  experimental  data  AT  10  F compared 
reasonably  with  the  theoretical  prediction.  However,  the  data  for 
AT  > 10 °F  are  far  below  the  theoretical  estimation.  The  reason  is  quite 
clear,  from  Fig.  12:  the  experiments  for  AT  > 10°F  did  not  follow  the 

slope  of  the  laminar  path  and  turned  toward  the  buoyancy-dominant  direc- 
tion. The  data  for  AT  > 10“F  are  affected  by  the  buoyancy  forces; 
therefore,  they  are  far  below  the  theoretical  predicted  values. 

t- 

The  study  of  the  developing  flow  in  the  region  of  0(aRe/Gr^)  is 
just  being  undertaken,  and  will  be  reported  separately  when  it  is  com- 
pleted. The  intuitive  expectation  on  the  flow  development  in  the  region 
of  0(aRe/Gr^)  would  be  that  the  unfavorable  pressure  gradient  developed 
on  the  top  wall  of  the  pipe  can  cause  flow  separation  and  possibly  trigger 
early  flow  transition.  Especially,  the  size  of  the  region  0(a)  shrinks 
when  the  heating  rate  increases. 

Personal  communication  from  Carl  Gazley,  Jr.,  The  Rand  Corporation. 
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Appendix  A 

NUMERICAL  INTEGRATION  OF  SINE  INTEGRALS 


The  integrals  of  Eqs.  (27)  and  (30a)  can  be  put  in  the  general  form  of 


f 


F(a)  sin  az  da 


(A-1) 


where  F(a)  may  have  integrable  singularity  at  a = 0.  The  integral  (A-1)  can 
be  separated  into  two  parts: 


From  0 - ir/z 


I 


sin  az 


where 


(A-2) 


G = limit  F(a) 
a 0 


(A-3) 


From  tt/z  - °° 


II 


sin  az  da 


1 

z 


00 


E 

k=l 


/" 

•\tt 


(k+1)  7T 


sin  a • da 


1. 

z 


sin  a 


F 


g + kir 
z 


da 


(A-4) 


The  infinite  series  in  (A-4)  is  summed  by  using  Euler's  transform, 


such  as 
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E (-1)’ 

k=l 


a + kir 


= E (-1)’ 

k=l 


a + kiT 


+ (-1) 


N+1 


M 

E v+1 

k=0  ^ 


(-1) 


A u. 


(A-5) 


where 


Uq  = F 


a + (N  + 1)  TT 


Au. 


k ^+1  \ 


(A-6) 


to  speed  the  convergence  of  the  summation.  N = 30  and  M = 10  are  used  to 
generate  data  which,  we  found,  are  far  more  than  sufficient. 

Equations  (A-2)  and  (A-3)  are  numerically  Integrated  by  using  Gaussian 
quadrature.  After  comparing  the  results,  we  found  that  either  ten  terms  or 
sixteen  terms  gives  sufficient  accuracy.  We  use  sixteen  terms. 

The  computer  program  is  listed  in  Appendix  B. 
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Appendix  B 
COMPUTER  PROGRAM 
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//Y1670LSy  JOB  ('^472, 150,  Y„  15)  „ • ISYB  0*  ,CIASS=D 
//COMP  EXEC  FOfiTCG,  P "GIOP.  GO=70K,LIBL='?YS1  .FOi-TEFF* 
//F05T.SYSIN  DD  * 

C*  PROGRAR  PIPE(FAST  FODEIER  TF&KSFOEM)  EOF,  5,  i97€ 

THE  ENTRANCE  FLOH  IN  A HEATED  PIPE 
CP  THE  SOLUTION  COVERS  Z=0(A) 


C 


c 


c 

c 

c 

c 

c 

c 

c 


(, 

c 

c 

L 

c 


IMPLICIT  HEAL*e  (A-li,0-Z)  , SEAL*  4 ($) 

DIMEHSIOH  G (50)  ,U  (10)  , H16  (16)  , X16  ( 16)  , HIO  (10|  , }f10  |10)  ,H  ( 16)  , XA  (16) 
COMMON  PAI 

DATA  NDGjNDU, NDH/50 , 10, 1 6/ 

DATA  H 16/0. 0271 5245941 1754D0, 

1 0.0  271  5245941  1754  1)0, 

2 0. 062253523938648DC, 

3 0.062253523938648D0, 

4 0.09515851  1682493DO, 

5 0.095158511682493D0, 

6 0. 124628971255534D0, 

7 0. 124628971255534D0, 

8 0. 149595988816577D0, 

9 0. 149595988816577D0, 

A 0. 169156519395003PO, 

B 0. 169156519395003D0, 

C 0. 182603415044924D0, 

D 0. 182603415044924EO, 

E 0. 1894506 1 045506BD0, 

F 0. 189450610455068D0/ 


DATA  X 16/0. 989400934991650P0 , 

1 -0.98940093499165010, 

2 0.944575023073^331)0, 

3 -0. 944575023073233L0, 

4 0. 865631202387S32DO, 

5 -0.86563120238783210, 

6 0. 755404408355003D0, 

7 -0. 75S404408355003PC, 

8 0.61787624440264410, 

9 -0. 61 7876244402644L0, 

A 0.458016777657227D0, 

B -0.458016777657227BC, 

C 0. 281603550779259D0, 

D -0. 281603550779259DO, 

E 0.095C  1 25098376  37t'0, 

F -0.09501250983763710/ 


DATA  Hi  0/0. 06  6 67  1 3 li  43  0 86  88  DO, 

1 0.0666^134430868880, 

2 0.14945134915058100, 

j 0.14845134915058100, 

4 0.21908636251598200, 

5 0 . 21  9086,  3625  15982D0  , 

6 0.26926671930999610, 

7 0. 26926671 9309996O0, 

8 0. 29552422471475300, 

9 0.  29  55242  2 47  1 475  3 1)0/ 


■ )A  :A  X 1 0/0  . 97  3006  52  85  1 7 1 72  I'O  , 

1 - 0 , 97  3906  5^  8“=  1 71  72  I'O  , 

2 0.865063^6  66889851)0, 

i -0.86504336668898500, 

4 0. 679409'^6R299024  I C„ 


0001 
00  02 

0007 
0004 
OOQp 
000-; 
000^ 

0008 

0009 

0010 
0011 
0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 
00  3 3 
00.34 
0035 
00  36 
00  37 
00  38 
0039 
004  0 
00  4 1 
0042 
00  4 3 

0044 

0045 

0046 

0047 
00  48 

0049 

0050 

0051 

0052 
00  5 3 
0054 
0 0 5- 
00  5 6, 

00^7 

0 0 5 fi 
0059 

or  ho 

0 0 6 ■) 
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C 5 -0. 67940  9568299  024  CO, 

C 6 0. 433395J94129247D0, 

C 7 -0. 433395394129247D0, 

C 8 0.  148874338981631  DO, 

C 9 -0.  148874338981531  DO/ 

C 

C*  M=NO.  OF  IEKH3  FOR  GAUSRIF.K  gUADSAOUPE 
C*  MG+NE=NO  OF  TERMS  FOR  THE  INFINITE  SERIES 

C»  H=  WSIGHIONG  FUNCTION  FOR  GAUSSIAN  QUADRATURE 

C»  XA=LOCATION  FOE  GAUSSIAN  QUADHATUPE 
CALL  EERSET  (208,256,-1,1,1) 

NC  = 0 

1001  NC=NC»1 

READ  1 01 , E,Z,NSE 
PRINT  204,  NC,R,Z 
PAI=3. 141592654D0 
SQHPAI  = DSQPT  (PAI) 

A=PAI/Z 

C 

C N=10 

C NG=20 

C NE  = 5 

C DO  10  1=1,  N 

C H(I)=H10(I) 

C 10  XA  (I)  =X1  0 (I) 

C GO  TO  1003 

1002  S=16 
NG  = 30 
NE  = 10 

DO  20  1 = 1, N 
H(I)  =H16  (I) 

20  XA  (I)  = X16  (I) 

NSS  = NSE*-1 
C 

C* 

C»  INTEGRATION  OF  P10 

C: 

1003  CALL  S IMP  (K,  A,F.,Z,P10  1 , 1 ,G,U,  H,  XA,  NDG,  NDU,  NDH) 

CALL  GAUSS (N, NG,NE, R, Z, PI  02,2 ,G,U, H, XA , NDG, NDU, NDH) 

P10  = 4. *Z#DSQBI  (A)  +P101+P1Q2 

P10=-1 ./SQEPAI*P10 

ASP10=-2. *DS0ET (2. *Z) » (1 . + 0. 2 5*B*?/4 . /Z**2) 

PRINT  201,  P101,P102 
C: 

C*  INTEGEATI01J  OF  P11  (CEOSS-FLOH  EFFECT) 

C* 

CALL  S IMP (F , A,R,Z, P 1 1 1 , 3,G,U, H, XA, NDG, NDU, NDH) 

CALL  GAUSS (N,NG,NE, S, Z, P1 1 2, 4,G,U, H, XA,NDG, NDU, NDH) 

PI  1 = R*Z*DSQRT (A) +P1 11  + P1 12 
PI 1 = -6 ./SQRPAI  + P1  1 

ASP1 1 = -3. *E* (D3QET ( 2 . *Z)  ♦ 0 . 5*  ( 1 . -0 . 2 5 *F» F) / (2 . *Z ) *♦ 1 . 5) 

PRINT  202,  P111,P112 
PRINT  203,  Pi  0, ASP1 0 , F 1 1 , ASP1 1 
101  FORMAT  (2D10. 3,110) 

204  FORMAT  (I10//5X,'(  R,  Z )=  ( ',D14.6,'  , ',D14.6,'  )') 

201  FORMAT  (/5 X , ' P 1 0 ( 0- A)  = ' , D 1 4 . 6 , 2 X , ’ Pi  0 ( A- X F)  = ' , D 1 4 . 6/) 

202  FORMAT  (5X,  ' PI  1 (0- A)  = ' , Dl  4 . t , 2 X , • P 1 ’ ( A- XF)  = ' , D 1 4 . 6/) 

203  FORMAT  (12X,****  RESULT  * ** • //2 X , ' DI SPL AC EM EN T EFFECT,  P10= 
1D14.6,'  ASP10=  ■ ,D14.6/4X , 'CEOSS-FLOH  EFFECT,  P11=  ', 
2D14.6,'  ASP11=  ',D14.6///) 

C 


0062 

0053 

0064 

0065 

006t 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 

0079 

0080 
008' 
0082 

0083 

0084 

0085 

0086 

0087 

0088 

0089 

0090 
009’ 

0092 

0093 

0094 

0095 

0096 

0097 

0098 
009° 
0100 
0101 
0102 

0103 

0104 

0105 

0106 
01  OP 
0108 
0109 
01  10 
01  n 
0112 
0113 
011“ 

0115 

0116 
01  17 
0118 
01  19 
0120 
0121 
01  22 
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c 

IF  ( NSE.EL;.1 

) GO  TO  1002 

0123 

c 

012U 

GO  TO  1001 

0125 

END 

0125 
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StlBaO'JTINE  SI»'P(N,A,E,7.,T0TAL,’JF,G,D,H,XA,NDG,KCU,NDH)  012-7 

c*  0128 

C*  NF  = 1 OE  3 0129 

I-1PLICIT  HEAL  + 8 (A-H.O-Z)  ,REAL*4  ($)  0130 

DIMENSION  G(NDG)  ,U  (NDU)  ,H  (KDH)  ,Xfl  (fiDH)  0131 

COMMON  PAI  0132 

A2=.5*A  0133 

TOTAL=0.  0134 

DO  10  1=1, N 0135 

Y=A2*(XA  (I) +1 .)  0136 

CALL  FUNC (Y,R,Z,FA, NF)  0137 

10  TOTAL=TOTAL+H (I) *FA  013« 

TOTAL=TOTAIFA2  0139 

EETUEN  0140 

EN  D 0141 


o o 
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SlJBHOU-INt  GAU5;S(N,NG,NF,,-,Z,Sa*n,>F,fJ,U,H,XJ,NDG,NDn,NDH)  0142 

C»  INTEGRATION  FROM  0.  TO  INFINITE  0143 

C*  N=NO.  OF  TERMS  FOE  GAUSS  0144 

C*  NG+NE=NO.  OF  TERMS  FO  3EEIES  0145 

C*  NF  = 3 OF  4 

IMPLICIT  EEAL*B  (A-H,0-Z)  ,F£AL*4  ($)  014"’ 

DIMENSION  G (NDG)  , LI  (NDO)  ,H  (NFH)  ,XA  (NDH)  0148 

COMMON  PAI  0149 

C***  GAUSSIAN  QUADRATURE  0150 

SUM1=0.  0151 

A2=.5*PAI  0152 

DO  50  1=1,  K 0153 

Y=A2*(1- +XA  (I) ) 0154 

G(I)=0.  0155 

DO  10  K=1,KG  0156 

¥Y= (Y+ K»PAI) /Z  0157 

CALL  FUNC (YY,P,Z,F, NF)  0159 

K1=MOD(K,2)  0159 

G(I)  =G  (I)  ♦ (-1  .)  **K1*F  0160 

GSUM=G(I)  0161 

IF  ( G(I).EQ.0.D0  ) GSUM=1.D0  0162 

EPR=DABS (F/GSUM)  0163 

IF  (ERE.  LE.  1 . D- 10)  GO  TO  50  0164 

10  CONTINUE  0165 

C*  EULER'S  TRANSFORM  0166 

DO  20  K=1,HE  0167 

YY=  (Y+ (NG«^  1 *K)  *PAI)  /Z  0168 

CALL  FUNC  (yY,R,2,F,HF)  0169 

20  U(K)=F  0170 

DO  30  L=2,NE  0171 

DO  30  K=L,NE  0172 

30  U (K)  =U  (K)  -U  (K- 1)  0173 

SUM=0.  0174 

DG=1.D0  0175 

DO  40  K=1,NE  0176 

KG=MOD (K, 2)  0177 

DG=2.*DG  0178 

40  SDH  = SUB- (-1  .)  *»KG/DG*U<K)  01''? 

C PRINT  201,  NF,I,F,2,G  (I)  ,SUM,ERR  0180 

KG  = MOD  (NG,  2) -s-l  0181 

G (I)  =G  <I)  + (-1  . ) "'*KG»SUM  0182 

50  SD«1=SUH1  4-H  (I)  *DSIN  (Y)  *G(I)  0183 

3UM1=SUM1 #A2/Z  0184 

RETURN  0185 

201  FORMAT  (/2X, 'GAUSS  MF , I , R , Z , G (I ) , S UM , EE E = •*'/  0186 

1 2I5,5D16.6/)  0187 

END  0188 


-51- 


SU3H0UTINE  FUNG  (XK,E,Z,F,N)  0189 

IMPLICIT  HtAL«8  (A-H,0-Z)  ,Hr,AL*4  ($)  0190 

RXK=H*XK  019^ 

ZXK=Z*XK  . 0192 

CALL  BESSEL  (XK , BO , B 1 , B 2)  0193 

CALL  BESSEL  (R XK , EBO , E t 1 , B B2)  0194 

X=XK  0195 

IF  (XK.EQ.O.)  X=1.  0196 

SXK=1. /DSOET  (X)  019"’ 

EXK=XK* (E-1 . DO)  0198 

FAC=DEXP <EXK)  0199 

IF  (N.GT.Z)  GO  TO  1003  0200 

F=FAC*BB0/B1  0201 

IF  (N.EQ.1)  F=F«DSIN (ZXK) -2.*Z  0202 

F=F*SXK  0203 

HE70RN  0204 

1003  F=FAC»BB1/ (B0+B2)  0205 

IF  (N.EO.3)  F=F*DSIN (ZXK) -0. 5»H*Z*XK**2  0206 

F=F#SXK/X**2  0207 

3ETUEK  0208 

END  0209 
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3UBR0UTINE  BESSEL  ( XK  , BC , B 1 , B2 ) 0210 

IMPLICIT  EEAL*R  (^-H,0-Z)  ,EEAL*4  ($)  0211 

C*  B(XK) =EESS£L  (XK)  »EXP (-XK)  0212 

IF  (XK . EQ. 0. ODO)  GO  TC  1002  0212 

KHAX=15  0214 

HHAX=10  0215 

B0=1.  021fc 

B1=1.  0217 

F0=1.  0218 

F1=1.  0219 

IF  (XK.GT.10.)  GO  TO  1001  0220 

X=XK*XK/4.  0221 

DO  10  K=1,KHAX  0222 

F0  = F0i‘X/K/K  0223 

F1=F1*X/K/(K+1)  0224 

B0=B0+F0  0225 

10  B1=B1+F1  0226 

B1=B1*XK»0. 5 0227 

B0  = B0*DEXP  (-XK)  0228 

B1  = B1*DEXP  (-XK)  0229 

GO  TO  30  0230 

1001  X=8.*XK  0231 

DO  20  K=2,NHAX  0232 

RK3=2.*K-3.  0233 

8K3=HK3**2  0234 

F0=F0*RK3/(K-1) /X  0235 

F1=F1*  (-1  . ) * (4.  -PK3)  / (K-1)  /X  02  36 

BO=BO+FO  0237 

20  B1=B1+F1  0238 

DEN=1. 00/DSQRT (2. *3. 14159*XK)  0239 

B0=B0*DEN  0240 

B1=B1*DEN  0241 

30  B2=B0-2./XK*B1  0242 

HSTORN  0243 

1002  B0=1.  0244 

B1=0.  0245 

B2=0.  0246 

RETURN  0247 

END  0248 

/♦  0249 

//GO.SYSIN  DD  f 0250 

.5  .5  0251 

.5  10.  0252 

/*  0253 
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1 


( R,  Z )= 
Pi  0 (0-A)  = 
Pi  1 (0-A)  = 

»♦: 

DISPLACEMENT 

CROSS-FLOH 


( 0.500000D 

-0.  1793 1 1D  01 
-0.316886D  00 
RESULT 

EFFECT,  P10= 
EFFECT,  P11= 


00  , 0.500000D  00  ) 

P10(A-XF)=  -0. 231 896D-01 
P11{A-XF)=  -0. 161 142D-03 

-0.180369D  01  ASP10=  -0.212500D  01 
-0.10U807D  01  ASP11=  -0.220313D  01 


2 


( B,  Z ) = ( 0. 500000D 

P10 (0-A) = -0. 566021D  01 

P11  (0-A)  = -0.712876D  00 

*♦*  RESULT  **» 

DISPLACEMENT  EFFECT,  P10= 
CROSS-FLOH  EFFECT,  P11= 


00  , 0.100000D  02  ) 

P10  (A-XF) = -0. 912196D 
P11  (A-XF) = -0.  109694D 

-0.894103D  01  ASP10= 
-0.670232D  01  ASP11= 


00 

00 

-0.894567D  01 
-0.671607D  01 


IHC900I  EXECUTION  TERMINATING  DUE  TO  EBROR  COUNT  FOE  ERROR  NUMBER  217 
IHC217I  FIOCS  - END  OF  DATA  SET  ON  UNIT  5 


TEACEBACK  ROUTINE  CALLED  FROM  ISN 

REG.  14 

REG.  15 

REG.  0 

REG.  1 

IflCOM 

0018C0F4 

001 8DF98 

00000003 

0018B9C8 

MAIN 

00014612 

0118B810 

FF00001 8 

001 9C758 

ENTRY  POINT=  0118B810 
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